/**
 * 
 */
package cn.cityhouse.avm.stat;

import org.apache.commons.math3.analysis.function.Sqrt;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;

import cn.cityhouse.avm.utils.MatrixUtil2;

/**
 * @author <a href="mailto:wangcl@lreis.ac.cn">ChenLiang Wang</a>
 *
 */
public class WeightedCovMatrices {

	public static RealMatrix estimateWeightedCov(RealMatrix xx, RealVector wt) {
		if (xx == null) {
			throw new IllegalArgumentException("null input matrix");
		}
		int n = xx.getRowDimension();
		if (wt != null) {
			if (wt.getDimension() != n) {
				throw new IllegalArgumentException("Length of wt must equal the number of row in x");
			}
			double sum_wt = 0;// =MatrixUtil2.getSum(wt);
			for (int i = 0; i < wt.getDimension(); i++) {
				double wti = wt.getEntry(i);
				if (wti < 0) {
					throw new IllegalArgumentException("weights must be non-negative and not all zero");
				}
				sum_wt += wti;
			}
			if (sum_wt == 0) {
				throw new IllegalArgumentException("weights must be non-negative and not all zero");
			}
			wt.mapDivideToSelf(sum_wt);
		} else {
			wt = MatrixUtil2.repeat(new ArrayRealVector(new double[] { 1 / xx.getRowDimension() }),
					xx.getRowDimension());
		}

		RealMatrix rlm_wtx = MatrixUtil2.productVectorByCol(xx, wt);
		RealVector rlv_center = MatrixUtil2.getColSum(rlm_wtx);
		RealMatrix rlm_cov = xx.copy();
		for (int i = 0; i < rlm_cov.getColumnDimension(); i++) {
			rlm_cov.setColumnVector(i, xx.getColumnVector(i).mapSubtract(rlv_center.getEntry(i)));
		}
		rlm_cov = MatrixUtil2.productVectorByCol(rlm_cov, wt.map(new Sqrt()));
		rlm_cov = rlm_cov.transpose().multiply(rlm_cov).scalarMultiply(1.0 / (1.0 - MatrixUtil2.getSumOfSquares(wt)));

		RealVector rlv_Is = MatrixUtil2.getReciprocal(MatrixUtil2.getDiagVector(rlm_cov).mapToSelf(new Sqrt()));
		RealMatrix rlm_rep = rlm_cov.copy();
		for (int c = 0; c < rlm_rep.getColumnDimension(); c++) {
			for (int r = 0; r < rlm_rep.getRowDimension(); r++) {
				rlm_rep.setEntry(r, c, rlv_Is.getEntry(c));
			}
		}
		RealMatrix rlm_cor = MatrixUtil2.getElementWiseProduct(MatrixUtil2.productVectorByCol(rlm_cov, rlv_Is),
				rlm_rep);
		return rlm_cor;
	}

	/**
	 * 
	 */
	public WeightedCovMatrices() {
		// TODO Auto-generated constructor stub

	}

}
